Airbnb_Project

Import packages

packages = c('sf', 'tidyverse', 'tmap', 'spatstat', 'raster', 'maptools', 'rgdal',
             'kableExtra', 'plotly', 'ggthemes', 'onemapsgapi', 'devtools')

# for each package, check if installed and if not, install it
for (p in packages){
  if(!require(p, character.only = T)){
    install.packages(p)
  }
  library(p,character.only = T)
}
Loading required package: sf
Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
Loading required package: tidyverse
Warning: package 'tidyverse' was built under R version 4.2.3
Warning: package 'ggplot2' was built under R version 4.2.3
Warning: package 'tibble' was built under R version 4.2.3
Warning: package 'tidyr' was built under R version 4.2.3
Warning: package 'readr' was built under R version 4.2.3
Warning: package 'purrr' was built under R version 4.2.3
Warning: package 'dplyr' was built under R version 4.2.3
Warning: package 'forcats' was built under R version 4.2.3
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.1     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.1     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Loading required package: tmap

Loading required package: spatstat

Loading required package: spatstat.data

Loading required package: spatstat.geom

spatstat.geom 3.0-5

Loading required package: spatstat.random

spatstat.random 3.1-3

Loading required package: spatstat.explore

Loading required package: nlme


Attaching package: 'nlme'


The following object is masked from 'package:dplyr':

    collapse


spatstat.explore 3.0-6

Loading required package: spatstat.model

Loading required package: rpart

spatstat.model 3.1-2

Loading required package: spatstat.linnet

spatstat.linnet 3.0-4


spatstat 3.0-3 
For an introduction to spatstat, type 'beginner' 


Loading required package: raster

Loading required package: sp


Attaching package: 'raster'


The following object is masked from 'package:nlme':

    getData


The following object is masked from 'package:dplyr':

    select


Loading required package: maptools

Checking rgeos availability: FALSE
Please note that 'maptools' will be retired during 2023,
plan transition at your earliest convenience;
some functionality will be moved to 'sp'.
    Note: when rgeos is not available, polygon geometry     computations in maptools depend on gpclib,
    which has a restricted licence. It is disabled by default;
    to enable gpclib, type gpclibPermit()

Loading required package: rgdal

Please note that rgdal will be retired during 2023,
plan transition to sf/stars/terra functions using GDAL and PROJ
at your earliest convenience.
See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
rgdal: version: 1.6-4, (SVN revision 1196)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
Path to GDAL shared files: C:/Users/65917/AppData/Local/R/win-library/4.2/rgdal/gdal
GDAL binary built with GEOS: TRUE 
Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
Path to PROJ shared files: C:/Users/65917/AppData/Local/R/win-library/4.2/rgdal/proj
PROJ CDN enabled: FALSE
Linking to sp version:1.6-0
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

Loading required package: kableExtra


Attaching package: 'kableExtra'


The following object is masked from 'package:dplyr':

    group_rows


Loading required package: plotly


Attaching package: 'plotly'


The following object is masked from 'package:raster':

    select


The following object is masked from 'package:ggplot2':

    last_plot


The following object is masked from 'package:stats':

    filter


The following object is masked from 'package:graphics':

    layout


Loading required package: ggthemes

Loading required package: onemapsgapi

Loading required package: devtools

Loading required package: usethis
pacman::p_load(sp, sf, rgdal, spNetwork, tmap)

Data

Geospatial Data

Import Data

beijing_sf <- st_read("data/Geospatial/beijing.geojson") %>%
  st_transform(crs=4555)
Reading layer `beijing' from data source 
  `C:\Quanfang777\IS415-GAA\Project\data\Geospatial\beijing.geojson' 
  using driver `GeoJSON'
Simple feature collection with 16 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 115.4172 ymin: 39.44247 xmax: 117.5071 ymax: 41.05861
Geodetic CRS:  WGS 84

Select Column for Analysis

beijing_sf <- select(beijing_sf, neighbourhood, geometry)
beijing_sf
Simple feature collection with 16 features and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 115.4172 ymin: 39.44247 xmax: 117.5071 ymax: 41.05861
Geodetic CRS:  New Beijing
First 10 features:
          neighbourhood                       geometry
1                东城区 MULTIPOLYGON (((116.4423 39...
2                西城区 MULTIPOLYGON (((116.3916 39...
3                昌平区 MULTIPOLYGON (((116.0427 40...
4       大兴区 / Daxing MULTIPOLYGON (((116.7347 39...
5                房山区 MULTIPOLYGON (((116.2466 39...
6      怀柔区 / Huairou MULTIPOLYGON (((116.279 40....
7  门头沟区 / Mentougou MULTIPOLYGON (((115.563 39....
8        密云县 / Miyun MULTIPOLYGON (((116.8826 40...
9       平谷区 / Pinggu MULTIPOLYGON (((117.3813 40...
10     延庆县 / Yanqing MULTIPOLYGON (((116.279 40....

Check CRS Again

st_crs(beijing_sf)
Coordinate Reference System:
  User input: EPSG:4555 
  wkt:
GEOGCRS["New Beijing",
    DATUM["New Beijing",
        ELLIPSOID["Krassowsky 1940",6378245,298.3,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Topographic mapping."],
        AREA["China - onshore."],
        BBOX[18.11,73.62,53.56,134.77]],
    ID["EPSG",4555]]

Check if there is invalid geometry

length(which(st_is_valid(beijing_sf ) == FALSE))
[1] 0

Check Missing Value

beijing_sf[rowSums(is.na(beijing_sf))!=0,]
Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  New Beijing
[1] neighbourhood geometry     
<0 rows> (or 0-length row.names)

Aspatial Data

beforecovid <- read_csv("data/aspatial/beijing_beforecovid.csv")
Warning: One or more parsing issues, call `problems()` on your data frame for details,
e.g.:
  dat <- vroom(...)
  problems(dat)
Rows: 25921 Columns: 106
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (51): listing_url, last_scraped, name, summary, space, description, expe...
dbl (39): id, scrape_id, host_id, host_listings_count, host_total_listings_c...
lgl (16): thumbnail_url, medium_url, xl_picture_url, host_is_superhost, host...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
aftercovid <- read_csv("data/aspatial/beijing_aftercovid.csv")
Rows: 6050 Columns: 18
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (5): name, host_name, neighbourhood, room_type, last_review
dbl (11): id, host_id, latitude, longitude, price, minimum_nights, number_of...
lgl  (2): neighbourhood_group, license

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Select Column for Analysis

beforecovid <- subset(beforecovid, select=c(id, street,neighbourhood_cleansed,latitude, longitude))
aftercovid <- subset(aftercovid, select=c(id,neighbourhood,latitude, longitude))

Check Missing Value

sum(is.na(beforecovid$latitude))
[1] 0
sum(is.na(aftercovid$latitude))
[1] 0
sum(is.na(beforecovid$longitude))
[1] 0
sum(is.na(aftercovid$longitude))
[1] 0

Visualize the listings

beforecovid_sf <- 

st_as_sf(beforecovid, 
                          # our coordinates are the longitude and latitude
                          coords = c("longitude", 
                                     "latitude"), 
                          # the geographical features are in longitude & latitude, in decimals
                          # as such, WGS84 is the most appropriate coordinates system
                          crs = 4326) %>%st_transform(crs=4555)
  #afterwards, we transform it to SVY21, our desired CRS
aftercovid_sf <- st_as_sf(aftercovid, 
                          # our coordinates are the longitude and latitude
                          coords = c("longitude", 
                                     "latitude"), 
                          # the geographical features are in longitude & latitude, in decimals
                          # as such, WGS84 is the most appropriate coordinates system
                          crs = 4326) %>% st_transform(crs=4555)

The listing before covid

tmap_mode("plot")+
  qtm(beijing_sf) +
  tm_shape(beforecovid_sf)+
  tm_dots()
tmap mode set to plotting

The listing after covid

tmap_mode("plot")+
  qtm(beijing_sf) +
  tm_shape(aftercovid_sf)+
  tm_dots()
tmap mode set to plotting

Network Analysis

RailWay

network <- st_read(dsn="data/Geospatial/Network/shape",
                   layer="railways")%>%
  st_transform(crs=4555)
Reading layer `railways' from data source 
  `C:\Quanfang777\IS415-GAA\Project\data\Geospatial\Network\shape' 
  using driver `ESRI Shapefile'
Simple feature collection with 8859 features and 3 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 116.0801 ymin: 39.68011 xmax: 116.77 ymax: 40.17997
Geodetic CRS:  WGS 84
network_sf <- select(network,osm_id,type, geometry)
network_sf
Simple feature collection with 8859 features and 2 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 116.0801 ymin: 39.68011 xmax: 116.77 ymax: 40.17997
Geodetic CRS:  New Beijing
First 10 features:
     osm_id    type                       geometry
1   9853311    rail LINESTRING (116.4297 39.900...
2   9853512    rail LINESTRING (116.4901 39.900...
3  23426905    rail LINESTRING (116.4913 39.900...
4  24818634 disused LINESTRING (116.1977 39.897...
5  24818659  subway LINESTRING (116.1485 39.947...
6  24818783    rail LINESTRING (116.4349 39.896...
7  24834552    rail LINESTRING (116.319 39.8942...
8  24850764    rail LINESTRING (116.4238 39.901...
9  24850768    rail LINESTRING (116.4241 39.900...
10 24851159    rail LINESTRING (116.4137 39.899...
plot(network)

st_crs(network_sf)
Coordinate Reference System:
  User input: EPSG:4555 
  wkt:
GEOGCRS["New Beijing",
    DATUM["New Beijing",
        ELLIPSOID["Krassowsky 1940",6378245,298.3,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Topographic mapping."],
        AREA["China - onshore."],
        BBOX[18.11,73.62,53.56,134.77]],
    ID["EPSG",4555]]
length(which(st_is_valid(network_sf) == FALSE))
[1] 0
network_sf[rowSums(is.na(network_sf))!=0,]
Simple feature collection with 0 features and 2 fields
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  New Beijing
[1] osm_id   type     geometry
<0 rows> (or 0-length row.names)

石景山区shijingshan

shijingshan_sf <- beijing_sf %>% filter(neighbourhood == "石景山区")
plot(shijingshan_sf)

network_shijingshan <- subset(network_sf, lengths(st_intersects(network_sf, shijingshan_sf))!=0,)
plot(network_shijingshan)

shijingshan_beforecovidlisting <- subset(beforecovid_sf, lengths(st_intersects(beforecovid_sf, shijingshan_sf))!=0,)
shijingshan_aftercovidlisting <- subset(aftercovid_sf, lengths(st_intersects(aftercovid_sf, shijingshan_sf))!=0,)

Before Covid Listing Visualization

tmap_mode("plot")
tmap mode set to plotting
tm_shape(shijingshan_sf) +
  tm_polygons() +
tm_shape(network_shijingshan) +
  tm_lines() +
tm_shape(shijingshan_beforecovidlisting) +
  tm_dots(size = 0.01,
          col = "blue",
          border.col="black",
          border.lwd=0.5)

After Covid Listing Visualization

tmap_mode("plot")
tmap mode set to plotting
tm_shape(shijingshan_sf) +
  tm_polygons() +
tm_shape(network_shijingshan) +
  tm_lines() +
tm_shape(shijingshan_aftercovidlisting) +
  tm_dots(size = 0.01,
          col = "blue",
          border.col="black",
          border.lwd=0.5)

Change them to spatial point, spatial line and spatial polygon

shijingshan_sf_as <- as_Spatial(shijingshan_sf)
shijingshan_beforecovidlisting_as <- as_Spatial(shijingshan_beforecovidlisting)
shijingshan_aftercovidlisting_as <- as_Spatial(shijingshan_aftercovidlisting)
network_shijingshan_as <- as_Spatial(network_shijingshan)
shijingshan_sp <- as(shijingshan_sf_as, "SpatialPolygons")
shijingshan_beforecovidlisting_sp <- as(shijingshan_beforecovidlisting_as, "SpatialPoints")
shijingshan_aftercovidlisting_sp <- as(shijingshan_aftercovidlisting_as, "SpatialPoints")
network_shijingshan_sl <- as(network_shijingshan_as, "SpatialLines")
shijingshan_sp <-spTransform(shijingshan_sp,
                        CRS("+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs"))
shijingshan_beforecovidlisting_sp <-spTransform(shijingshan_beforecovidlisting_sp,
                        CRS("+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs"))
shijingshan_aftercovidlisting_sp <-spTransform(shijingshan_aftercovidlisting_sp,
                        CRS("+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs"))
network_shijingshan_sl <-spTransform(network_shijingshan_sl,
                        CRS("+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs"))

Network Constrained KDE (NetKDE) Analysis

Preparing the lixels objects

network_shijingshan_sf1  <- st_as_sf(network_shijingshan_sl )
crs <- "+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs"
network_shijingshan_sf1 <- st_transform(network_shijingshan_sf1, crs)
lixels <- lixelize_lines(network_shijingshan_sf1, 
                         700, 
                         mindist = 350)
shijingshan_beforecovidlisting_sf1 <- st_as_sf(shijingshan_beforecovidlisting_sp)
shijingshan_aftercovidlisting_sf1 <- st_as_sf(shijingshan_aftercovidlisting_sp)
crs <- "+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs"
shijingshan_beforecovidlisting_sf1<- st_transform(shijingshan_beforecovidlisting_sf1, crs)
shijingshan_aftercovidlisting_sf1<- st_transform(shijingshan_aftercovidlisting_sf1, crs)

Generating line centre points

samples <- lines_center(lixels)

Performing NetKDE for Beijing Before Covid

densities <- nkde(network_shijingshan_sf1, 
                  events = shijingshan_beforecovidlisting_sf1,
                  w = rep(1,nrow(shijingshan_beforecovidlisting_sf1)),
                  samples = samples,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 100,
                  sparse = TRUE,
                  verbose = FALSE)
samples$density <- densities
lixels$density <- densities
samples$density <- samples$density*100
lixels$density <- lixels$density*100
tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(lixels)+
  tm_lines(col="density")+
tm_shape(shijingshan_beforecovidlisting_sf1)+
  tm_dots()

Network Constrained G- and K-Function Analysis

# kfun_shijingshan_bef<- kfunctions(network_shijingshan_sf1, 
#                              shijingshan_beforecovidlisting_sf1,
#                              start = 0, 
#                              end = 500, 
#                              step = 50, 
#                              width = 50, 
#                              nsim = 20, 
#                              resolution = 25,
#                              verbose = FALSE, 
#                              conf_int = 0.05)

Performing NetKDE for Beijing After Covid

densities1 <- nkde(network_shijingshan_sf1, 
                  events = shijingshan_aftercovidlisting_sf1,
                  w = rep(1,nrow(shijingshan_aftercovidlisting_sf1)),
                  samples = samples,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 100,
                  sparse = TRUE,
                  verbose = FALSE)
samples$density1 <- densities
lixels$density1 <- densities
samples$density1 <- samples$density*100
lixels$density1 <- lixels$density*100
tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(lixels)+
  tm_lines(col="density1")+
tm_shape(shijingshan_aftercovidlisting_sf1)+
  tm_dots()
# kfun_shijingshan_aft<- kfunctions(network_shijingshan_sf1,
#                              shijingshan_aftercovidlisting_sf1,
#                              start = 0,
#                              end = 500,
#                              step = 50,
#                              width = 50,
#                              nsim = 20,
#                              resolution = 25,
#                              verbose = FALSE,
#                              conf_int = 0.05)

Haidian 海淀

haidian_sf <- beijing_sf %>% filter(neighbourhood == "海淀区")%>%
  st_transform(crs=4555)
plot(haidian_sf)

network_haidian <- subset(network_sf, lengths(st_intersects(network_sf, haidian_sf  ))!=0,)%>%
  st_transform(crs=4555)
plot(network_haidian )

haidian_beforecovidlisting <- subset(beforecovid_sf, lengths(st_intersects(beforecovid_sf, haidian_sf))!=0,) %>%
  st_transform(crs=4555)
haidian_aftercovidlisting <- subset(aftercovid_sf, lengths(st_intersects(aftercovid_sf, haidian_sf))!=0,)%>%
  st_transform(crs=4555)

Before Covid Listing Visualization

tmap_mode("plot")
tmap mode set to plotting
tm_shape(haidian_sf) +
  tm_polygons() +
tm_shape(network_haidian) +
  tm_lines() +
tm_shape(haidian_beforecovidlisting) +
  tm_dots(size = 0.01,
          col = "blue",
          border.col="black",
          border.lwd=0.5)

After Covid Listing Visualization

tmap_mode("plot")
tmap mode set to plotting
tm_shape(haidian_sf) +
  tm_polygons() +
tm_shape(network_haidian) +
  tm_lines() +
tm_shape(haidian_aftercovidlisting) +
  tm_dots(size = 0.01,
          col = "blue",
          border.col="black",
          border.lwd=0.5)

crs <- "+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs"
network_haidian<- st_transform(network_haidian, crs)
haidian_beforecovidlisting<- st_transform(haidian_beforecovidlisting, crs)
haidian_aftercovidlisting<- st_transform(haidian_aftercovidlisting, crs)

Network Constrained KDE (NetKDE) Analysis

lixels_haidian <- lixelize_lines(network_haidian, 
                         700, 
                         mindist = 350)
samples_haidian <- lines_center(lixels_haidian)
samples_haidian<- st_transform(samples_haidian, crs)
densities_haidian <- nkde(network_haidian, 
                  events = haidian_beforecovidlisting,
                  w = rep(1,nrow(haidian_beforecovidlisting)),
                  samples = samples_haidian,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 10,
                  sparse = TRUE,
                  verbose = FALSE)
samples_haidian$density <- densities_haidian
lixels_haidian$density <- densities_haidian
samples_haidian$density <- samples_haidian$density*100
lixels_haidian$density <- lixels_haidian$density*100
tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(lixels_haidian)+
  tm_lines(col="density")+
tm_shape(haidian_beforecovidlisting)+
  tm_dots()

After Covid Visualization

densities_haidian <- nkde(network_haidian, 
                  events = haidian_aftercovidlisting,
                  w = rep(1,nrow(haidian_aftercovidlisting)),
                  samples = samples_haidian,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 10,
                  sparse = TRUE,
                  verbose = FALSE)
samples_haidian$density <- densities_haidian
lixels_haidian$density <- densities_haidian
samples_haidian$density <- samples_haidian$density*100
lixels_haidian$density <- lixels_haidian$density*100
tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(lixels_haidian)+
  tm_lines(col="density")+
tm_shape(haidian_aftercovidlisting)+
  tm_dots()

FENGTAI 丰台

fengtai_sf <- beijing_sf %>% filter(neighbourhood == "丰台区 / Fengtai")%>%
  st_transform(crs=4555)
plot(fengtai_sf)

network_fengtai <- subset(network_sf, lengths(st_intersects(network_sf, fengtai_sf  ))!=0,)%>%
  st_transform(crs=4555)
plot(network_fengtai)

fengtai_beforecovidlisting <- subset(beforecovid_sf, lengths(st_intersects(beforecovid_sf, fengtai_sf))!=0,) %>%
  st_transform(crs=4555)
fengtai_aftercovidlisting <- subset(aftercovid_sf, lengths(st_intersects(aftercovid_sf, fengtai_sf))!=0,)%>%
  st_transform(crs=4555)

Before Covid Visualization

tmap_mode("plot")
tmap mode set to plotting
tm_shape(fengtai_sf) +
  tm_polygons() +
tm_shape(network_fengtai) +
  tm_lines() +
tm_shape(fengtai_beforecovidlisting) +
  tm_dots(size = 0.01,
          col = "blue",
          border.col="black",
          border.lwd=0.5)

After Covid Visualization

tmap_mode("plot")
tmap mode set to plotting
tm_shape(fengtai_sf) +
  tm_polygons() +
tm_shape(network_fengtai) +
  tm_lines() +
tm_shape(fengtai_aftercovidlisting) +
  tm_dots(size = 0.01,
          col = "blue",
          border.col="black",
          border.lwd=0.5)

Network Constrained KDE (NetKDE) Analysis

Preparing the lixels objects

crs <- "+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs"
network_fengtai<- st_transform(network_fengtai, crs)
fengtai_beforecovidlisting<- st_transform(fengtai_beforecovidlisting, crs)
fengtai_aftercovidlisting<- st_transform(fengtai_aftercovidlisting, crs)
lixels_fengtai <- lixelize_lines(network_fengtai, 
                         5600, 
                         mindist = 2800)
samples_fengtai <- lines_center(lixels_fengtai)
densities_fengtai <- nkde(network_fengtai, 
                  events = fengtai_beforecovidlisting,
                  w = rep(1,nrow(fengtai_beforecovidlisting)),
                  samples = samples_fengtai,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 10,
                  sparse = TRUE,
                  verbose = FALSE)
samples_fengtai$density <- densities_fengtai
lixels_fengtai$density <- densities_fengtai
samples_fengtai$density <- samples_fengtai$density*100
lixels_fengtai$density <- lixels_fengtai$density*100
tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(lixels_fengtai)+
  tm_lines(col="density")+
tm_shape(fengtai_beforecovidlisting)+
  tm_dots()
# kfun_childcare <- kfunctions(network_fengtai, 
#                              fengtai_beforecovidlisting,
#                              start = 0, 
#                              end = 1000, 
#                              step = 50, 
#                              width = 50, 
#                              nsim = 50, 
#                              resolution = 50,
#                              verbose = FALSE, 
#                              conf_int = 0.05)

After covid

densities_fengtai <- nkde(network_fengtai, 
                  events = fengtai_aftercovidlisting,
                  w = rep(1,nrow(fengtai_aftercovidlisting)),
                  samples = samples_fengtai,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 10,
                  sparse = TRUE,
                  verbose = FALSE)
samples_fengtai$density <- densities_fengtai
lixels_fengtai$density <- densities_fengtai
samples_fengtai$density <- samples_fengtai$density*100
lixels_fengtai$density <- lixels_fengtai$density*100
tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(lixels_fengtai)+
  tm_lines(col="density")+
tm_shape(fengtai_aftercovidlisting)+
  tm_dots()